Bayesian Statistics

A Bayesian is one who, vaguely expecting a horse, and catching a glimpse of a donkey, strongly believes he has seen a mule

Why might someone use Bayesian rather than traditional frequentist?

• If they have prior beliefs as to the truth before an experiment is conducted
• If you want to predict how the distribution of a population might change, rather than a point estimate. For example, traditionally, we might try and predict the mean number of people who choose to order after going on our website.

Quick recap of the ‘Classical Approach’ (‘Frequentist’)

set.seed(0)
p <- 0.5
n <- 10

Let’s say we have a coin, and we want to know the probability that we flip heads. Frequentists assume that there is only one true, fixed, unknown probability that the coin yields heads. For ease of notation, let’s call the constant parameter \(p(\text{head}) = \theta\).

To try and discover the true probability of getting a head, we flip the coin 10 times, and count the number of heads to estimate the probability. Frequentists believe that if we repeated this sampling proceedure of 10 flips many, many times, we should expect the average across all of our samples to be reflective of the true population value (hence why it is called ‘frequentist’).

# Create repeated samples
num_repeated_samples <- 10^3
repeated_samples <- rbinom(n = num_repeated_samples, size = n, prob = p)
repeated_samples.cum_p_mean <- cumsum(repeated_samples)/cumsum(rep(n, num_repeated_samples))

# Cumulative mean tends to p
plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = 1:length(repeated_samples), y = repeated_samples.cum_p_mean, name = 'Cumulative mean with sampling') %>%
  layout(showlegend = F
         , yaxis = list(range = c(0, 1))
         , title = "Cumulative sample means tends to true p") # , 

k <- round(n*p*0.5,0) # make up fake k that is 'unlikely'
sample_mean <- k/n

binomial_probability <- function(k, n, p) {
  n_choose_k <- factorial(n)/(factorial(k)*factorial(n-k))
  # NOTE - equivalent to gamma_function(k+n-k)/(gamma_function(k)*gamma_function(n-k), or 1/B(k,n-k)
  prob_sample <- n_choose_k*(p)^k*(1-p)^(n-k)
  return(prob_sample)
}

k_heads_maxprob <- round(binomial_probability(k,n,k/n),4)*100
k_heads_actprob <- round(binomial_probability(k,n,p),4)*100
np_heads_actprob <- round(binomial_probability(n*p,n,p),4)*100

But in the real world, it’s not possible to take lots of samples over and over again - and we generally just take one. Taking a sample one time, we might observe \(2\) heads. This is equivalent to creating a sample of size \(10\) and observing a mean probability of \(0.2\).

In frequentist statistics, we assume there is only one true probability for flipping heads (\(\theta\)). We thus want to find that probability that is most likely to give us the result we observed in our sample (a ‘maximum likelihood estimate’). On the left graph below, for every possible value that the probability of flipping heads could be on the x axis (somewhere between 0 and 1), the likelihood of observing \(2\) heads is plotted on the y axis, giving the orange line. For example, the likelihood of observing \(2\) heads if the probablity of flipping heads is \(0.2\) is \(30.2\%\). If the probablity of flipping heads is actually \(0.5\), then the likelihood of flipping 2 heads in our sample is only \(4.39\%\). It appears then, based solely on our sample, that our true probability is most likely to be \(0.2\).

We say estimate, however, because despite treating the true \(\theta\) as a fixed (non-random) value, frequentists acknowledge there is random sampling error. Even if the probability of flipping heads remaining constant, sometimes there will be 2 heads, sometimes there will be 10 heads etc in our sample, simply due to chance. To clarify, frequentists believe the probability is fixed (at \(\theta\)), but the number of heads is a random variable conditioned by that fixed parameter: if \(\theta\) is genuinely \(0.2\), then observing \(2\) heads is more likely than observing \(10\) heads, but both results are possible.

We can illustrate how the number of heads observed follows a probability distribution conditioned by the fixed value of \(\theta\) by plotting the likelihood of observing different numbers of heads. In the graph on the right, the probability of getting different numbers of heads between 0 and 10 is plotted in green if the value of theta is \(0.2\). (A bayesian might also consider other values of theta, such as the blue one if the probability is \(0.5\), but will discuss that later).

# for all the different values of theta, p(H=2|θ)
theta_prob_k_heads <- sapply(seq(0,1,0.01), FUN = function(p) binomial_probability(k=k,n=n,p=p))
# for all the different values of outcome k, p(H=k|θ=0.2)
likelihood_est <- sapply(seq(0,10,1), FUN = function(s) {binomial_probability(k=s,n=n,p=k/n)})
# for all the different values of outcome k, p(H=k|θ=0.5)
likelihood_act <- sapply(seq(0,10,1), FUN = function(s) {binomial_probability(k=s,n=n,p=p)})


subplot(
  plot_ly(type = 'scatter', mode = 'lines') %>%
    add_trace(x = seq(0,1,0.01), y = theta_prob_k_heads, name = paste0('P(H=',k,'|θ)')) %>%
    layout(xaxis = list(title = 'θ'))
  , plot_ly(type = 'scatter', mode = 'lines+markers') %>%
    add_trace(x = seq(0,10,1), y = likelihood_est, name = paste0('p(H=k|θ=',k/n,')')) %>%
    add_trace(x = seq(0,10,1), y = likelihood_act, name = paste0('p(H=k|θ=',p,')'), visible = 'legendonly') %>%
    layout(xaxis = list(title = '# heads'))
  , shareY = T
  , nrows = 1)

# 
# plot_ly(type = 'scatter', mode = 'lines') %>%
#   add_trace(x = x, y = likelihood_est, name = 'Estimated from sample') %>%
#   add_trace(x = x, y = likelihood_act, name = 'Actual likelihood') %>%
#   layout(title = "Likelihood function: X ~ Binom(n,p)"
#          , xaxis = list(title = 'theta')
#          , yaxis = list(title = 'likelihood'))

Thus, to acknowledge that there is randomness associated with this sample mean, frequentists construct confidence intervals. Again to follow the frequentist interpretation, a 95% confidence interval means that, if the sampling proceedure was repeated many, many times, that 95% of the time the true value of \(\theta\) would fall within the confidence interval. (Nuanced note - this is not the same as saying that, in one sample, there is a 95% chance that the true value of \(\theta\) lies within it. This is because \(\theta\) is a fixed value, not a random variable, so it either falls within the sample confidence interval or not.)

true.p.se <- sqrt(p*(1-p)/n)

## BINOMIAL (NORMAL APPROXIMATION)
# The larger the number of successes, np, and failures, n(1-p), the better the normal approximation for confidence intervals (p near 0.5, n as large as possible)
sample.binom.k <- rbinom(n = num_repeated_samples, size = n, prob = p)
sample.binom.p <- sample.binom.k/n
  ## population standard error/confidence interval
    true.binom.p.conf.02.5 <- sample.binom.p-1.96*true.p.se
    true.binom.p.conf.97.5 <- sample.binom.p+1.96*true.p.se
    true.binom.theta_in_conf <- true.binom.p.conf.02.5 < p & true.binom.p.conf.97.5 > p
  ## sample standard error/confidence interval
    sample.binom.p.se <- sqrt(sample.binom.p*(1-sample.binom.p)/n)
    sample.binom.p.conf.02.5 <- sample.binom.p-1.96*sample.binom.p.se
    sample.binom.p.conf.97.5 <- sample.binom.p+1.96*sample.binom.p.se
    sample.binom.theta_in_conf <- sample.binom.p.conf.02.5 < p & sample.binom.p.conf.97.5 > p

## BETA DISTRIBUTION (CONTINUOUS VERSION OF BINOMIAL)
# Since the binomial distribution is non-continuous, we can get confidence intervals closer to 95% from the equivalent beta distribution for the same np & n(1-p)
sample.beta.k <- rbeta(n = num_repeated_samples, shape1 = p*n, shape2 = (1-p)*n)*n
sample.beta.p <- sample.beta.k/n
  ## population standard error/confidence interval
    true.beta.p.conf.02.5 <- sample.beta.p-1.96*true.p.se
    true.beta.p.conf.97.5 <- sample.beta.p+1.96*true.p.se
    true.beta.theta_in_conf <- true.beta.p.conf.02.5 < p & true.beta.p.conf.97.5 > p
  ## sample standard error/confidence interval
    sample.beta.p.se <- sqrt(sample.beta.p*(1-sample.beta.p)/n)
    sample.beta.p.conf.02.5 <- sample.beta.p-1.96*sample.beta.p.se
    sample.beta.p.conf.97.5 <- sample.beta.p+1.96*sample.beta.p.se
    sample.beta.theta_in_conf <- sample.beta.p.conf.02.5 < p & sample.beta.p.conf.97.5 > p

# First one-time sample confidence intervals
one_time.sample.p.se <- sqrt((k/n)*(1-k/n)/n)
one_time_sample.conf_02.5 <- (k/n)-1.96*one_time.sample.p.se
one_time_sample.conf_97.5 <- (k/n)+1.96*one_time.sample.p.se

# One time sample either contains true value of theta or not
plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = rep(c(one_time_sample.conf_02.5, one_time_sample.conf_97.5), each = 3), y = c(1,-1,0,0,-1,1), name = 'Confidence interval') %>%
  add_trace(x = p, y = 0, mode = 'markers', name = 'True value', marker = list(size=50, symbol = 'x')) %>%
  layout(showlegend = F
         , yaxis = list(range = c(-2, 2), zeroline = F, showticklabels = F)
         , xaxis = list(range = c(one_time_sample.conf_02.5 - 0.1, one_time_sample.conf_97.5 + 0.1))
         , title = "One-time sample confidence interval either contains true value of theta or not")

# Cumulative confidence intervals contain 95% of the time
plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = 1:num_repeated_samples
            , y = cumsum(true.beta.theta_in_conf)/cumsum(rep(1,num_repeated_samples))
            , name = 'Cumulative mean with sampling') %>%
  layout(showlegend = F
         #, yaxis = list(range = c(0, 1))
         , title = "Cumulative confidence intervals tends to contain p 95% of the time")

NA

In order to obtain these results, certain conditions have to be met, including:

• The randomly sampled observations are independent, so observing any one value does not have an impact on what subsequent values we observe might be. In our example, flipping a head once does not impact the likelihood of flipping it again. (If running an experiment in practice, this might mean that our population is sufficiently large when compared to our sample that we don’t need to be concerned with replacement. Another thing to watch out for is network effects - imagine an airline increases price of seats as they sell. Say we want to experiment in offering some customers a special offer for airline tickets, our variant, and some not, our control: the customers with the special offer might be more likely to buy, and this will shorten the supply, increasing the price for those without the special offer, and decreasing their likelihood to buy).
• Any sampled observations are identically distributed to the population (e.g. if we take our sample from the USA, they might have a different propensity to watch TV than the UK.)
• Because it is a random sample, the values observed are characterised by a probability distribution associated with \(\theta\). (e.g. we don’t know what the number of heads oberved will be before sampling, but the likelihood of what we observe is pre-determined through its associated probability distribution with \(\theta\))

Okay, so what exactly does a Bayesian approach mean?

Bayesian inference is broadly similar to frequentist inference. The key difference though is that, under a Bayesian approach, the true value of \(\theta\) is treated as a random variable, rather than a fixed constant.

So, for example, we might believe that in our economy, \(90\)% of coins are fair, and \(10\)% of coins are biased, with biased coins flipping heads \(20\)% of the time.

Under the guise of being frequentists previously, we had no prior belief as to what \(\theta\) would be: we only made an educated guess based on the 2 heads we observed. We believed that there is one fixed value of the true \(\theta\), and that there was only randomness associated with sampling, so we created a confidence interval to portray this.

Now, as bayesians, we are also incorporating our prior beliefs as to what \(\theta\) is - that it is either fair or biased. We believe that there is not just random error associated in sampling - we also treat \(\theta\) as an additional random variable, since we don’t know whether our coin is fair or biased either, and fair coins will have a different sample mean distribution than biased ones.

fair_coin_k_prob <- binomial_probability(k, n, p)
biased_coin_k_prob <- binomial_probability(k, n, k/n)

Given that we observe 2 heads, the likelihood of this occuring if the coin is fair is \(4%\)%, and \(30%\)% if the coin is biased. In the absence of informative prior information then, it might appear that the coin is more likely to be biased.

However, we also need to take into account that \(90\)% of coins in the population are fair, so even before observing 2 heads, we have some belief about whether the coin is fair or biased (i.e. the probability distribution associated with the random variable \(\theta\)).

fair_coin_prob <- (fair_coin_k_prob * 0.9)/((fair_coin_k_prob * 0.9) + (biased_coin_k_prob * 0.1))
biased_coin_prob <- (biased_coin_k_prob * 0.1)/((fair_coin_k_prob * 0.9) + (biased_coin_k_prob * 0.1))

# data.frame(
# matrix(data = c(binomial_probability(2,10,0.5)*0.9 # p(fair, 2 heads)
#                 , (1-binomial_probability(2,10,0.5))*0.9 # p(fair, not 2 heads)
#                 , 0.9 # p(fair)
#                 , binomial_probability(2,10,0.2)*0.1 # biased, 2 heads
#                 , (1-binomial_probability(2,10,0.2))*0.1 # biased, not 2 heads
#                 , 0.1
#                 , binomial_probability(2,10,0.5)*0.9 + binomial_probability(2,10,0.2)*0.1
#                 , (1-binomial_probability(2,10,0.5))*0.9 + (1-binomial_probability(2,10,0.2))*0.1
#                 , 1
#                 )
#        , nrow = 3
#        , dimnames = list(c('Flip 2 heads', 'Not flip 2 heads', 'Total')
#                          , c('Fair', 'Biased', 'Total'))
#        )
# )

..and so we can find the likelihood we got a fair coin, given the fact we observed 2 heads:

\[ P(\text{fair}|H=2) =\frac{P(\text{H=2 and fair})}{P(\text{H=2 and fair})+P(\text{H=2 and biased})} =\frac{P(\text{H=2|fair}) \times P(\text{fair})}{P(\text{H=2})} \]

\[ P(\text{fair}|H=2) =\frac{0.0439 \times 0.9}{0.0439 + 0.302}=0.567 \]

\[ P(\text{biased}|H=2) =\frac{P(\text{H=2 and biased})}{P(\text{H=2 and fair})+P(\text{H=2 and biased})} =\frac{P(\text{H=2|biased}) \times P(\text{biased})}{P(\text{H=2})} \] \[ P(\text{biased}|H=2) =\frac{0.302 \times 0.1}{0.0439 + 0.302}=0.433 \]

In other words, our prior belief (before sampling and observing any data) was that the likelihood of our coin having \(\theta\) equal to \(0.5\) was \(90\%\), and the likelihood of it having theta equal to \(0.2\) was \(10\%\). After observing 2 heads, we update our beliefs, so that the likelihood of the probability of our coin flipping heads (\(\theta\)) being equal to \(0.5\) is now equal to \(56.7\%\), and of it being equal to \(0.2\) being \(43.3\%\)

It is this incorporation of prior data that is either seen as Bayesian’s biggest advantage or pitfall. On the one hand, experiments are not abstract devices, and some knowledge about the process being investigated before obtaining the data is known and arguably should be incorporated. On the other, incorporating subjective opinions, particularly strong ones, may mean that you do not learn the true values you are trying to derive. Bayesian is thus analysis that uses a set of observations to change opinion rather than as a means to determine ultimate truth.

To help describe bayesian analysis, the following terms are often used:

• Prior distribution: \(Pr(\theta)\) - represents existing belief about \(\theta\) (Represents what was thought before seeing the data). For example, in the binomial coin example above, we have prior knowledge that 90% of coins in our economy are fair, so \(P(\theta)=P(Biased)=0.1\)
• Evidence: \(X\) - what we just observed (\(2\) heads)
• Marginal probability: \(Pr(X)\) - the total probability of the data across all possible values of the parameter \(\theta\). It doesn’t actually depend on \(\theta\) and isoften referred to as the proportionality factor/normalising constant (it makes sure all the scenarios are modelled. For example, \(P(H=2) = P(H=2|Fair) \times P(Fair) + P(H=2|Biased) \times P(Biased)\))
• Likelihood function: \(Pr(X|\theta)\) - the probability of \(\theta\) given the data observed (Represents the new data available). For example, \(P(H=2|Biased)=0.3019899\)
• Joint probability density function: \(Pr(X,\theta)=Pr((X|\theta).Pr(\theta)\) - used to modify prior beliefs through Bayes Theorem
• Posterior distribution: \(Pr(\theta|X)\) - the posterior density represents the knowledge about the model parameters after observing the data (Represents what is now thought given both prior data and data just obtained)
• Conjugancy: occurs when the posterior distribution is in the same family of probability density functions as the prior belief, but with new parameter values (updated to reflect what has been learned from the data). The posterior comes from the same family as the prior, rather than the data’s distribution.

And thus, we can then describe Bayes Theorem:

\[ P(\theta|X) =\frac{P(X|\theta) \times P(\theta)}{P(X|\theta) \times P(\theta) + P(X|\theta') \times P(\theta')} =\frac{P(X|\theta) \times P(\theta)}{P(X)} \]

Thinking about Bayesian in terms of predicting getting a head next time:

There are two sources of uncertainty when building statistical models to predict things:

• Uncertainty due to the fact any future value is itself a random event, \(P(y|\theta)\) (because future events are essentially samples that have error conditioned by \(\theta\))
• Uncertainty in the parameter values which have been estimated on the basis of past data, \(P(\theta|X)\) (for example, a coefficient estimated in a regression model)

Classical models deal with point 1, making a point estimate in the future based on the ‘best’ parameters estimated from past data. This is equivalent to forming a likelihood (i.e. incorporating new data). What they don’t do though is think about point 2, that the seemingly optimum model estimated itself might be incorrect since the true values of \(\theta\) are not fixed, and this is equivalent to not building in any prior beliefs. In a way, by giving a point estimate, the models generate a false sense of precision, and the confidence intervals arounds any estimate is only predicated on the idea there is sampling randomness, rather than randomness associated with \(\theta\) itself.

In other words, classical models seek to best fit the new evidence (traditionally through maximum likelihood estimation for generalized linear models), whereas bayesian models incorporate prior beliefs as well.

To predict the likelihood of getting heads next time, given our evidence, we do the following:

\[ \begin{aligned} P(\text{heads}|k=2) \\ = P(\text{heads}|\text{fair},k=2) \times P(\text{fair}|k=2) + P(\text{heads}|\text{biased},k=2) \times P(\text{biased}|k=2) \\ = 0.5 \times P(\text{fair}|k=2) + 0.2 \times P(\text{biased}|k=2) \\ = 0.5 \times \frac{P(k=2|\text{fair}) \times P(\text{fair})}{P(k=2)} + 0.2 \times \frac{P(k=2|\text{biased}) \times P(\text{biased})}{P(k=2)} \\ = 0.5 \times 0.567 + 0.2 \times 0.433 \\ = 0.3701 \end{aligned} \]

Note that we have subbed in our results for the probability that we have a fair or a biased coin using bayes theorem, to allow us to predict the probability we get heads next time.

If we chose a purely frequentist approach, we would have guessed that the likelihood is \(0.2\), since we observed \(2\) heads, and can only take one value (and this was the maximum likelihood estimate). If we chose to only use our prior informtion that 90% of coins were fair, we would have guessed \(0.9 \times 0.5 + 0.1 \times 0.2 = 0.47\). However, since we have incorporated both the prior information and the new evidence, we predict a result that is somewhere in the middle, at \(0.3701\).

We can formalise this:

\[ f(y|x)= \int f(y|\theta,x)f(\theta|x) \partial(\theta) \] In other words, for all our potential estimates of what \(\theta\) could be (\(\delta(\theta)\)), how likely each of those values for \(\theta\) might be given the number of heads we have observed (\(f(X|\theta)\)), and how likely we were to observe our number of heads given those parameters \(P(y|\theta,x)\), we can make a prediction for what our number of heads might be in the next flip.

This now allows us to tackle problems where the prior comes from continuous distributions.

Continuous prior beliefs

In the last example, we had a very strong prior - \(\theta\) could only take 2 values, and was 90% likely to be fair, so it had a very large influence on the result. However, Bayesian tends to most valuable where both the new data and past prior have more equal influence, and the posterior forming a distribution that sits somewhere between that estimated by the prior and likelihood function.

So - let’s say the last time we tried flipping 10 times, we got heads 5 times (our prior) and the new coin flipped heads 2 times (our new evidence). We can then get the following equations:

Our prior is a binomial:

\[ \text{Prior: } P(\theta) \sim \text{Binom}[n_{1},k_{1}/n_1] = \text{Binom}[10,0.5] \] Our likelihood is also binomial:

\[ \text{Likelhood: } P(X|\theta) \sim\text{Binom}[n_{2},k_{2}/n_2] = \text{Binom}[10,0.2] \] Our posterior thus follows a binomial distribution too:

\[ \text{Posterior: } P(\theta|X) = \frac{P(X|\theta)P(\theta)}{\int P(X|\theta)P(\theta) d\theta} \approx \text{Binom}\left[n_1+n_2,\frac{k_1+k_2}{n_1+n_2}\right] = \text{Binom}[7,0.35] \\ \]

In other words, we obtain our posterior estimate by multiplying our likelihood function and prior beliefs, and dividing it by all possible probabilities given our evidence that theta could take. In practice, this denominator just ensures that our probabilities sum to one.

We can plot this to see what is looks like:

# prior <- sapply(x, FUN = function(theta) {beta_probability(theta,k1+1,n1-k1+1)/(n1+1)})
# likelihood_f <- sapply(x, FUN = function(theta) {beta_probability(theta,k2+1,n2-k2+1)/(n2+1)})
# posterior <- sapply(x, FUN = function(theta) {beta_probability(theta,k1+k2+2,n1-k1+1+n2-k2+1)/(n1+n2+2)})

# prior <- dbeta(x, k1+1, n1-k1+1)/(n1+1)
# likelihood_f <- dbeta(x, k2+1, n2-k2+1)/(n2+1)
# posterior <- dbeta(x, k1+k2+1, n1-k1+n2-k2+1)/(n1+n2+1)

prior <- dbinom(seq(0,10,1), n1, k1/n1)
likelihood_f <- dbinom(seq(0,10,1), n2, k2/n2)
posterior <- dbinom(seq(0,20,2), n1+n2, (k1+k2)/(n1+n2))*2

plot_ly(type = 'scatter', mode = 'lines+markers') %>%
  add_trace(x = seq(0,10,1), y = prior, name = 'Prior belief: Horse') %>%
  add_trace(x = seq(0,10,1), y = likelihood_f, name = 'Likelihood (New data): Donkey') %>%
  add_trace(x = seq(0,10,1), y = posterior, name = 'Posterior Belief: Mule') %>%
  layout(title = "Updating Prior Beliefs")

The posterior curve is both taller and narrow than both the distributions estimated from the prior and the likelihood functions. This is because, by including both prior data and new evidence, our sample size has increased, so this reflects greater confidence that the observed number of heads,conditioned by the probability distribution of \(\theta\), will lie within a smaller interval.

However, in general, we might have a continuous prior or likelihood function, and the beta distribution is a continuous functional form of the binomial distribution.

A quick introduction to the beta distribution: it can model lots of different functional forms using \(a\) and \(b\) as parameters, by inputting them into the following function:

\[ \text{Beta}[a,b] \sim \theta^{a-1}\times(1-\theta)^{b-1} \]

Now let’s define the binomial distribution for comparison:

\[ \text{Binom}(n,p) = {n \choose k } [\theta^{k}\times(1-\theta)^{n-k}] \]

Now to get geeky - you might have already noticed how eerily similar the beta distribution is to the binomial distribution from the formula - and in fact you can use the beta distribution as a continous function equivalent to the binomial function.

\[ p(\text{H=}k| \theta) \sim \theta^k\times(1-\theta)^{n-k} \approx \text{Beta}[k+1,n-k+1] \]

Both binomial and beta likelihoods are scaled though (normalized), so that the sum of their probabilites total 1. For the beta distribution, the scaling is achieved like this:

\[ Beta[\alpha,\beta] = \theta^{\alpha-1}\times(1-\theta)^{\beta-1} \times \frac{1}{\text{B}(\alpha,\beta)} \] Where

\[ \frac{1}{\text{B}(\alpha,\beta)} = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \times \Gamma(\beta)} = \frac{(\alpha+\beta-1)!}{(\alpha-1)! \times (\beta-1)!} \]

So if successes \(\alpha = k+1\) and failures \(\beta = n - k + 1\), then:

\[ \frac{(\alpha+\beta-1)!}{(\alpha-1)! \times (\beta-1)!}=\frac{((k+1)+(n-k+1)-1)!}{((k+1)-1)! \times ((n-k+1)-1)!}=\frac{(n+1)!}{(k)!\times(n-k)!} = {n+1 \choose k } = (n+1) \times {n \choose k } \]

And thus:

\[ Beta[k+1,n-k+1] = (n+1) \times \text{Binom}(n,p) \]

This relationship helps us combine binomial likelihoods and beta priors to create a posterior in what is called conjugancy. Conjugancy occurs when the posterior distribution is in the same family of probability density functions as the prior belief, but with new parameter values (updated to reflect what has been learned from the data). The posterior comes from the same family as the prior, rather than the data’s distribution.

So - if we want to use a beta distribution for our binomial prior, we can model this using a beta distribution by setting \(a\) and \(b\) to both equal to 1

\[ \text{Prior: } p(\theta) \approx \text{Beta}[k_{1}+1,n_{1}-k_{1}+1] = \text{Beta}[6,10-5+1] = \text{Beta}[5,6] \]

\[ \text{Likelihood: } p(k|\theta)=\text{Binom}[n_{2},k_{2}/n_2] = \text{Beta}[k_{2}+1,n_{2}-k_{2}+1] = \text{Beta}[2+1,10-2+1] = \text{Beta}[3,9] \]

\[ \text{Posterior: } p(\theta|k) \approx \\ \text{Beta}[5,6] \times \text{Beta}[3,9] = \\ \text{Beta}[k_1+k_{2}+2,n_1-k_1+n_{2}-k_{2}+2] = \text{Beta}[9,15] \]

# prior <- sapply(x, FUN = function(theta) {beta_probability(theta,k1+1,n1-k1+1)/(n1+1)})
# likelihood_f <- sapply(x, FUN = function(theta) {beta_probability(theta,k2+1,n2-k2+1)/(n2+1)})
# posterior <- sapply(x, FUN = function(theta) {beta_probability(theta,k1+k2+2,n1-k1+1+n2-k2+1)/(n1+n2+2)})

# prior <- dbeta(x, k1+1, n1-k1+1)/(n1+1)
# likelihood_f <- dbeta(x, k2+1, n2-k2+1)/(n2+1)
# posterior <- dbeta(x, k1+k2+1, n1-k1+n2-k2+1)/(n1+n2+1)

x <- seq(0, 1, length=100+1)

prior <- dbeta(x, k1+1, n1-k1+1)#/(n1+1)
likelihood_f <- dbeta(x, k2+1, n2-k2+1)#/(n2+1)
posterior <- dbeta(x, k1+k2+2, n1-k1+n2-k2+2)#/((n1+1)*(n2+1))
# prior <- prior/sum(prior)
# likelihood_f <- likelihood_f/sum(likelihood_f)
# posterior <- posterior/sum(posterior)

plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = x, y = prior, name = 'Prior belief: Horse') %>%
  add_trace(x = x, y = likelihood_f, name = 'Likelihood (New data): Donkey') %>%
  add_trace(x = x, y = posterior, name = 'Posterior Belief: Mule') %>%
  layout(title = "Updating Prior Beliefs")

Why might we want to use beta though rather than just binomial as our prior? Well a beta distribution is really flexible, so we can use it to model a range of different priors depending on previous experiments, or even as a way of regularising through the bayesian paradigm.

So - if we want to create a uniform distribution for our prior, we can model this using a beta distribution by setting \(a\) and \(b\) to both equal to 1

\[ p(\theta) \sim \text{Beta}[1,1] \sim \theta^{1-1}\times(1-\theta)^{1-1}=\theta^{0}\times(1-\theta)^{0}= 1 \text{ (for all values of }\theta) \]

And we can model our likelihood using the binomial function

\[ p(X| \theta) = p(\text{H=}k| \theta) = \text{Binom}(n,p) \sim \theta^k\times(1-\theta)^{n-k} \]

And finally - the posterior is proportional to the prior multiplied by the likelihood:

\[ p(\theta) \times p(X| \theta) \sim [\theta^{a-1}\times(1-\theta)^{b-1}] \times [\theta^k\times(1-\theta)^{n-k}] = \theta^{a-1+k}\times(1-\theta)^{b-1+n-k} \]

Which we can purely parametize in a beta distribution (you can see that, the stronger the prior - and hence the larger a and b - the smaller the impact new evidence k and n has on the posterior distribution):

\[ \theta^{a-1+k}\times(1-\theta)^{b-1+n-k} \sim \text{Beta}[a+k,b+n-k] \]

gamma_function <- function(n) {
  return(factorial(n-1))
}

beta_probability <- function(theta,a,b) {
  numerator <- ((theta)^(a-1))*((1-theta)^(b-1))
  denominator <- (gamma_function(a)*gamma_function(b))/gamma_function(a+b)
  return(numerator/denominator)
}

a <- 1
b <- 1

# uniform distribution
# prior <- dbeta(x = x, shape1 = 1, shape2 = 1)
x <- seq(0, 1, length=100+1)
prior <- sapply(x, FUN = function(theta) {beta_probability(theta,a,b)})

# binomial distribution
likelihood_f.binom <- sapply(x, FUN = function(theta) {binomial_probability(k,n,theta)})
likelihood_f.binom.int <- likelihood_f.binom
likelihood_f.binom.int[round(x*n,0) != x*n] <- NA

# beta distribution
likelihood_f.beta <- sapply(x, FUN = function(theta) {beta_probability(theta,k+1,n-k+1)*1/(n+1)})

posterior <- sapply(x, FUN = function(theta) {beta_probability(theta,a+k+1,b+n-k+1)*1/(n+1)})

plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = x, y = prior/sum(prior), name = 'Prior: Uniform = Beta(1,1)') %>%
  add_trace(x = x, y = likelihood_f.binom.int, name = 'Likelihood: ~ Binom(k=2,p=theta,n)', mode = 'markers') %>%
  add_trace(x = x, y = likelihood_f.beta, name = 'Likelihood: ~ Beta(p=theta,k,n-k)') %>%
  add_trace(x = x, y = posterior, name = 'Posterior: ~ Beta(a+k,b+n-k)') %>%
  layout(xaxis = list(title = 'Theta'))

Ignoring 91 observationsIgnoring 91 observations

As you can see, with a very ‘objective’ prior, the Bayesian posterior is very similar to the result derived from the frequentist evidence-only distribution (the Bayesian likelihood function).

Markov-chain Monte-Carlo (MCMC): Estimating posterior distributions when the prior/likelihood density functions are non-conjugant

Sometimes though, we have priors and likelihood functions that do not play nicely with each other, and hence we cannot define the posterior density function analytically. In this instance, we can make many draws from the posterior distribution for various values of \(\theta\) to infer what the posterior distribution looks like.

This is the ‘monte-carlo’ bit - we randomly choose values of \(\theta\) to make draws from the posterior from.

However, the aim of this is to create a histogram that most closely resembles the posterior distibution that we are trying to estimate. For this histogram to be reflective of the posterior distribution, we need to take draws more frequently from the areas higher density from the posterior, and less frequently where they are of lower density.

This is where the ‘markov-chain’ bit comes in - when making a new draw, whether we add this new value to the histogram, or the old one again to the histogram, depends on the probability density associated with the new point vs the old one.

The instructions for the MCMC algorithm (using the Metropolis-Hastings method) is as follows:

\[ \text{for(i in num_draws) } \{\\ \hspace{3cm}\theta_{new,i}:= X \sim N(\theta_{recorded,i-1},\sigma^2) \\ \hspace{3cm}posterior_{new,i}:= P(X|\theta_{new,i}) \times P(\theta_{new,i}) \\ \hspace{3cm}posterior_{old,i}:= P(X|\theta_{old,i}) \times P(\theta_{old,i}) \\ \hspace{3cm}r_{num,i}:= r_i \sim \text{unif}(0,1)\\ \hspace{3cm} \text{if}\left(\frac{posterior_{new,i}}{posterior_{old,i}} > r_{num,i}\right) \{\\ \hspace{6cm}\theta_{recorded,i}:=\theta_{new,i}\\ \hspace{6cm}\theta_{recorded,i}:=\theta_{old,i}\\ \hspace{3cm}\} \\ \} \]

What does this mean?

• With each value of \(\theta_{recorded,i-1}\) we record from the previous iteration, we take a random walk from that point to trial another candidate \(\theta_{new,i}\) (drawing from a normal distribution with mean given as the previously recorded \(\theta_{recorded,i-1}\), and fixed standard deviation \(\sigma^2\))
• If that area is of higher density than the previous one (so the ratio is 1 or greater) then we record that
• However, we don’t automatically reject it. If it is bigger than a number randomly generated from a uniform distribution (\(r_i\)) then we still accept it (to allow if to explore the distribution rather than get stuck in local maxima)
• We can then make a histogram of all the recorded \(\theta_{recorded}\) across all iterations, and it should mirror the true distribution of the posterior.

Let’s take our example before, where our prior had 5 of 10 flips, and our evidence was 2 of 10 flips, to check how well this works. We run the algorithm 10,000 times (in general use, we might exclude a lot first draws, called the burn, but in this instance as long as we take enough samples it should still create a reflective histogram).

As a reminder:

• The prior was \(\text{Beta}[k1+1,n1-k1+1]=\text{Beta}[6,6]\)
• The likelihood was \(\text{Binom}[2,0.2]\)
• So we expect the posterior \(\sim \text{Beta}[k1+k2+2,n1+n2-k1-k2+2]=\text{Beta}[8,15]\)):

## MCMC M-H
theta_v <- c()
for(i in 1:10^5) {
  if(i == 1) {
    theta_old <- rnorm(1,p,0.05)
    theta_new <- rnorm(1,p,0.05)
  } else {
    theta_old <- chosen_theta
  }
  theta_new <- rnorm(1,theta_old,0.05)
  if(theta_new < 0 | theta_new > 1) {
    break
  }
  posterior_new <- dbeta(theta_new,k1+1,n1-k1+1) * dbinom(k2,n2,theta_new) # new prior * likelihood
  posterior_old <- dbeta(theta_old,k1+1,n1-k1+1) * dbinom(k2,n2,theta_old) # old prior * likelihood
  if(posterior_new/posterior_old > runif(1)) {
    chosen_theta <- theta_new
  } else {
    chosen_theta <- theta_old
  }
  theta_v <- c(theta_v,chosen_theta)
}

plot_ly() %>%
  add_trace(x=theta_v,type='histogram',name='frequency') %>%
  add_trace(x=density(theta_v)$x, y=density(theta_v)$y,type='scatter',mode='lines',yaxis='y2',name='Density approximation') %>%
  add_trace(x=seq(0,1,0.01), y=dbeta(seq(0,1,0.01),k1+k2+2,n1+n2-k1-k2+2),type='scatter',mode='lines',yaxis='y2',name=paste0('True posterior: Beta(',k1+k2+2,',',n1-k1+n2-k2+2,')')) %>%
  layout(yaxis2=list(overlaying="y",side="right",rangemode = "tozero"))

Thinking about Bayesian in terms of predicting getting a head next time (continuous prior):

So now we are set up to make a prediction for getting a heads next time, given our prior beliefs and our new evidence.

Let’s remind ourselves of the equation to make predictions:

\[ f(y|x)= \int f(y|\theta,x)f(\theta|x) \partial(\theta) \]

We now know \(p(\theta|x)\) - our posterior distribution. So now we just need to estimate \(p(y|\theta,x)\). So how do we do this?

Well we know that \(y\) follows a bernoulli distribution, and we can solve this through maximum likelihood estimation. In particular, for binary predictions, we often assume that the likelihood function follows a logistic distribution, given as \(p(y|\theta,X)=(1+e^{-X\beta})^{-1}\) - You can see this in another note I have written here.

Great, so just need to integrate the product of \(f(y|\theta,x)\) and \(f(\theta|x)\) with respect to \(\theta\).

y_dist <- c()
theta_options <- seq(0,1,0.01)
for(i in theta_options) {
  theta_posterior_dist = dbeta(i,k1+k2+2,n1-k1+n2-k2+2)
  y_logit = (1+exp(-(k2)*theta_posterior_dist))^(-1)
  y_dist <- c(y_dist,theta_posterior_dist*y_logit)
}

plot_ly() %>%
  add_trace(x=seq(0,1,0.01), y=y_dist,type='scatter',mode='lines',name=paste0('P(y|X)')) %>%
  layout(title='Likelihood of observing heads next time')

NA

This is more interesting if we had varying explanatory variables \(X\) - and then we can derive different posterior distributions for \(y\) given different values for those explanatory features.

Bayesian Credible intervals and HPD (need to change this to \(y\))

So now we have our posterior distribution for \(\y\), how might we summarise our new understanding as to where \(\y\) might lie in the range of?

This leads on to Bayesian Credible intervals - the Bayesian equivalent to frequentist’s confidence intervals. Because Bayesian estimates distributions rather than fixed values, credible intervals are arguably more intuitive than confidence intervals: they reflect the a (95%) probability that of here the value generated by the random variable theta will fall within the range (rather than the frequentist interpretation that the true value will fall within the confidence interval 95% of the time).

Like a confidence interval, a credible interval is symmetric around the mean, and spans the x-axis enough to cover 95% of the area under the probability distribution curve. And that is good if the posterior distribution is also symmetric - but as can be seen in our example above, that is not necessarily the case, where the posterior is skewed.

If the distribution is skewed then, it could be better to compute the region of highest posterior density (HPD) - which by definition will be the narrowest interval across the potential values for theta, but cover 95% of the area under the probability curve. In other words, we still maintain the same 95% probability of theta lying within the region, but the region is the narrowest it possibly can be.

So to do this, we want to find the interval between \(\theta_{low}\) and \(\theta_{high}\) such at:

\[ \int_{\theta_{low}}^{\theta_{high}} p(\theta|x) d\theta= 1-\alpha \]

But finding this can be tricky, particularly with computing closed intervals. The code below shows a brute force solution. The way to think about this is as follows:

• Credible intervals: Given the mean of the distribution, expand the range outwards equally until 50% of the area under the curve is captured.
• Highest posterior density: Imagine a horizontal line moving slowly towards the x axis, intercepting the distribution twice. The line stops descending once it is low enough that the auc is 50%.

plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = x, y = y_dist, name = 'Beta(9,15)') %>%
  add_trace(x = rep(0.36,2), y = c(0,max(sapply(highest_posterior_density, FUN = function(x) {dbeta(x,k1+k2+2,n1-k1+n2-k2+2)}))), name = 'Mean') %>%
  add_trace(x = pmin(pmax(x, credible_intervals[1]),credible_intervals[2]), y = likelihood_f, name = 'Credible Interval', fill = 'tozeroy', line = list(dash = 'dot')) %>%
  layout(title = 'Credible Interval (equidistant around mean)')


plot_ly(type = 'scatter', mode = 'lines') %>%
  add_trace(x = seq(0,1,0.001), y = sapply(seq(0,1,0.001), FUN = function(x) {dbeta(x,k1+k2+2,n1-k1+n2-k2+2)}), name = 'Beta(9,15)') %>%
  add_trace(x = rep(0.36,2), y = c(0,max(sapply(highest_posterior_density, FUN = function(x) {dbeta(x,k1+k2+2,n1-k1+n2-k2+2)}))), name = 'Mean') %>%
  add_trace(x = highest_posterior_density, y = sapply(highest_posterior_density, FUN = function(x) {dbeta(x,k1+k2+2,n1-k1+n2-k2+2)}), name = 'Highest Posterior Density', fill = 'tozeroy', line = list(dash = 'dot')) %>%
  layout(title = 'Highest Posterior Density (smallest spread that contains 95% of data)')

The width of the credible interval is larger, ranging from 0.1883926 to 0.5616074 (range of 0.3732147), whereas width the range of the HPD is smaller, ranging from 0.186 to 0.56 (range of 0.374)

Some Bayesians argue though that credible intervals and HPDs are only really something that are useful to compare to frequentist interpretations, as the distribution is already given by the posterior.